線形回帰を実装してみよう

データ分析

#機械学習

じょんすみす

2017.12.01

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

こんにちは、小澤です。

当エントリは「Machine Learning Advent Calendar 2017」の1日目のエントリです。

今回は、初回ということでウォーミングアップがてら、線形回帰を実装してみましょう。

前提知識

さて、最初に前提となる知識について、解説していきます。

機械学習で扱うデータ

機械学習ではデータに含まれる性質を明らかにします。そのため、データがなければ機械学習はできません。

では、どのような形式でどのようなデータを用意すればいいのでしょうか？これは以下のようなテーブル形式になります。

sepal_length	sepal_width	petal_length	petal_width	species
5.1	3.5	1.4	0.2	Iris-setosa
4.9	3.0	1.4	0.2	Iris-setosa
4.7	3.2	1.3	0.2	Iris-setosa
4.6	3.1	1.5	0.2	Iris-setosa

これは、iris(あやめ)に関する観測データになっており、以下のような項目が含まれています。

sepal_length : がくの長さ
sepal_width : がくの幅
petal_length : 花びらの長さ
petal_width : 花びらの幅
species : 品種

このデータは機械学習における入門的なデータになっており、4つの要素から、あやめの品種を予測するようなモデルを作成する練習問題として利用されます。

データの中には、数値以外にも性別や都道府県などといった特定のどれかに当てはまるものや文字列で表現されるものもあります。機械学習を行う際は、基本的に対象を数値として扱う必要があるため、それらを何らかの形で変換してやる必要があります。

例えば、都道府県であれば、対応する都道府県コードを振ってやったり、以下のように各都道府県を列方向に展開して元データで対応する部分は「1」、それ以外は「0」を入れてやるなどで対応します。

変換前

名前住んでる都道府県

A 北海道

B 北海道

C 青森

D 秋田

E 秋田

... ...

F 沖縄
変換後

名前北海道青森秋田 ... 沖縄

A 1 0 0 ... 0

B 1 0 0 ... 0

C 0 1 0 ... 0

D 0 0 1 ... 0

E 0 0 1 ... 0

... ... ... ... ... ...

F 0 0 0 ... 1

名前	住んでる都道府県
A	北海道
B	北海道
C	青森
D	秋田
E	秋田
...	...
F	沖縄

名前	北海道	青森	秋田	...	沖縄
A	1	0	0	...	0
B	1	0	0	...	0
C	0	1	0	...	0
D	0	0	1	...	0
E	0	0	1	...	0
...	...	...	...	...	...
F	0	0	0	...	1

この形式は、One-Hotなベクトルなどと呼ばれます。都道府県コードなどは便宜的に数値を当てはめているだけですので、この値に対して四則演算などを行ってしまうと困るような手法を使う際に活用します。

また、文字列なども一度単語に分割したのち、元データの各文章における単語の出現回数に対応した上記のような形式に変換したりすることもあります(これはBag of Words; BoWなどと呼ばれたりもします)。

機械学習の種類

続いて、機械学習の種類を見ていきましょう。機械学習は、やりたいことに応じていくつかの種類に分けられます。

機械学習は大きく分けて以下のように分類できます。

教師あり学習
- 回帰
- 分類
教師なし学習
- クラスタリング, 次元削減など
その他
- 強化学習など

回帰

回帰は、様々な要因から推測される数値を予測するものになります。

例えば、web広告, 雑誌広告, テレビCMなどに出した数とその月の商品の売り上げが上記のようなテーブル形式で与えられているとします。このとき、過去の売り上げ実績を元に各広告の出稿数を変化させた時に売り上げがどのように変わるかをみるといった場面で利用できます。

分類

分類は、数値ではなくいくつかの項目のうちどれに当てはまるかを予測するようなものになります。例えば、最初に例として出したirisのデータの場合、4つの数値からどの品種であるかを予測するなどになります。

分類は、webアクセス者の年齢, 性別, 居住地域などから広告がクリックされるか否かのどちらかを当てるといった2つのうちどちらかとなるものはニ値分類、3つ以上のどれかのうちどれに当てはまるかは多値分類と呼ばれます。

教師なし学習

教師なし学習は、分類や回帰と異なり、データ中に正解となる値が存在しないようなものを対象とした手法になります。

例えば、スーパーに買い物に来たお客さんが実際に購入した商品の傾向を元にいくつかのグループに分けてみる、といった場面で利用できます。グループの別れ方は機械側に任せっきりになるので、人間にとって解釈するのが難しかったり期待したものとは異なる性質で分かれる可能性もあります。

その他

その他の手法に関してはここでは割愛します。

線形回帰とは

さて、今回の実装対象である線形回帰とはどういったものなのでしょうか？回帰というからには、数値を予測するものになります。

線形回帰は

[latex] y = ax + b [/latex]

のように観測されたデータxに対して、予測対象となる数値yが直線の関係で表されるモデルを作成するものになります。どのような計算を行うのかについてのより詳細な去年Alteryx Advent Calendarで解説したものがありますので、そちらをご覧ください。

線形回帰をしてみる #alteryx #07 | Alteryx Advent Calendar 2016 ｜ Developers.IO

実装してみる

予測してみる数値

今回は、Boston house-pricesを使ってみたいと思います。これはボストン市内の住宅の価格を予測するような問題になります。

データセット中の各列に関する解説は以下をご参照ください

scikit-learn に付属しているデータセット

まずはこのデータを読み込んで確認してみます。

from sklearn.datasets import load_boston
import pandas as pd

boston = load_boston()
boston_df = pd.DataFrame(boston.data, columns=boston.feature_names)
boston_df['target'] = boston.target

以下のようなデータを取得できます。

なお、ここではデータの取得にscikit-learn, 表示にPandasを利用していますが、以降の実装ではnumpyのみを利用していきます。

実装する

では、順に実装していきましょう。

データの前処理

さて、最初に切片の項(x_0)を1としてしまうことで、別枠で計算しなくて良くなります。そのため、切片付きのデータに変換しましょう。

#　現在のboston.dataよりも1列多いゼロ行列を作成
boston_data = np.zeros([boston.data.shape[0], boston.data.shape[1] + 1])

# 各要素にboston.dataの値を入れていく
for i in range(boston.data.shape[0]):
    for j in range(boston.data.shape[1]):
        boston_data[i, j] = boston.data[i, j]
    # 追加した切片項には1を入れる
    boston_data[i, boston.data.shape[1]] = 1

続いて、データの正規化を行います。データの正規化は各列の平均、標準偏差を使って値を変換します。

[latex] z = \frac{x - \mu}{\sigma} [/latex]

各列の値をこの変換後の値で更新していきます。

# 各列の平均
boston_data_mu = [sum(boston_data[:, row]) / len(boston_data[:, row]) for row in range(boston_data.shape[1])]

# 各列の標準偏差
boston_data_sigma = [math.sqrt(sum((boston_data[:, row] - boston_data_mu[row]) ** 2  / len(boston_data[:, row]))) for row in range(boston_data.shape[1])]

# 正規化
for i in range(boston_data.shape[0]):
    # 切片の項は平均1, 標準偏差0となるが、この計算をする際に0割のエラーが出るので、計算対象から除く
    for j in range(boston_data.shape[1]-1):
        boston_data[i, j] = (boston_data[i, j] - boston_data_mu[j]) / boston_data_sigma[j]

最後に初期値としてランダムなweightの値を設定します。

# 元データに含まれる列名リストのサイズに切片項を足した数分のw
w = random(len(boston.feature_names) + 1)

学習の処理

さて、いよいよ学習の処理を実装します。今回は最急勾配方で計算をしています。

alpha = 0.01
for i in range(50):
    y_hat = np.dot(boston_data, w)
    # 損失関数 : error = sum(((boston.target - y_hat) ** 2) / 2)
    w = w + alpha * (np.dot((boston.target - y_hat) , boston_data) / boston_data.shape[0])

まず最初にイテレーションとなるので、ループ内で処理を行います。損失関数の値が一定以下としてもいいのですが、今回は簡単のため、固定回数にています。

次に現在のweightの値で予測値を求めています。

y_hat = np.dot(boston_data, w)

各行の予測値は

[latex] \hat{y} = {\bf w}^T {\bf x} [/latex]

で求められます。これを全データでまとめて処理をするために、行列X(今回の場合boston_data)にたいして

[latex] {\bf \hat{y}} = X {\bf w} [/latex]

の計算を行っています。これの損失関数の勾配を計算するのですが、求めるのが勾配だけであれば損失関数そのもは計算しなくてもいいので、微分した値のみを計算します。

損失関数は

[latex] E = \frac{1}{2m} \sum_{i}^{m}(y_i - {\bf w}^T {\bf x}_i)^2 [/latex]

微分すると、

[latex] \frac{\partial E}{\partial {\bf w}_n} = - \frac{1}{m} \sum_{i}^{m} (y_i - {\bf w}^T {\bf x}_i) {\bf x}_{in} [/latex]

となります。

-np.dot((boston.target - y_hat) , boston_data) / boston_data.shape[0]

では、各列nごとに計算するのではなく、以下のように行列演算でまとめて計算しています。

[latex] - \frac{1}{m} (y - \hat{y})^T X [/latex]

この値で、wを更新しています。

[latex] {\bf w}_n = {\bf w}_n - \alpha \frac{\partial E}{\partial {\bf w}_n} \\ = {\bf w}_n + \frac{\alpha}{m} \sum_{i}^{m} (y_i - {\bf w}^T {\bf x}_i) {\bf x}_{in} [/latex]

w = w + alpha * (np.dot((boston.target - y_hat) , boston_data) / boston_data.shape[0])

結果を確認する

さて、最後に学習結果を確認してみます。誤差の値や決定係数など調べてもいいのですが、今回は見た目でどうなってるかわかる感じにしてみました。

学習のループ処理を以下のように変更しています。

for i in range(50):
    y_hat = np.dot(boston_data, w)
    # 10イテレーションごとにグラフを出力
    if (i % 10 == 0):
        result = pd.DataFrame({'y' : boston.target, 'y_hat' : y_hat})
        print(result.plot.scatter('y', 'y_hat'))
    # 損失関数 : error = sum(((boston.target - y_hat) ** 2) / 2)
    w = w + alpha * (np.dot((boston.target - y_hat) , boston_data) / boston_data.shape[0])

横軸に正解の値、縦軸に予測値をプロットしています。以下のように、学習が進むごとにいい感じの結果(だんだん正の相関がありそうな感じ)になっていってることがわかります。

学習したweightの値

最後に学習したweightの値を見てみましょう。

print(w)

[-0.62901299  0.65535304 -0.26547212  0.88373226 -0.34828324  2.1784461
 -0.31073857  0.5024018  -0.46448136 -0.15024068 -0.77601093  0.79928798
 -1.35060284  9.49582607]

weightがプラスになっているものは予測値を上げる効果を持ち、マイナス担っているものは下げる効果になっています。また、各値の大きさがどれだけ要因として大きいものなのかも確認できます。

今回は、この値のみですが、実際のライブラリを利用すれば、本当は影響していない(weightがゼロでも問題ない)要素がわかったりもします。